Network Pharmacology and Molecular Docking Elucidate the Underlying Pharmacological Mechanisms of the Herb Houttuynia cordata in Treating Pneumonia Caused by SARS-CoV-2

Used in Asian countries, including China, Japan, and Thailand, Houttuynia cordata Thumb (H. cordata; Saururaceae, HC) is a traditional herbal medicine that possesses favorable antiviral properties. As a potent folk therapy used to treat pulmonary infections, further research is required to fully elucidate the mechanisms of its pharmacological activities and explore its therapeutic potential for treating pneumonia caused by SARS-CoV-2. This study explores the pharmacological mechanism of HC on pneumonia using a network pharmacological approach combined with reprocessing expression profiling by high-throughput sequencing to demonstrate the therapeutic mechanisms of HC for treating pneumonia at a systemic level. The integration of these analyses suggested that target factors are involved in four signaling pathways, including PI3K-Akt, Jak-STAT, MAPK, and NF-kB. Molecular docking and molecular dynamics simulation were applied to verify these results, indicating a stable combination between four metabolites (Afzelin, Apigenin, Kaempferol, Quercetin) and six targets (DPP4, ELANE, HSP90AA1, IL6, MAPK1, SERPINE1). These natural metabolites have also been reported to bind with ACE2 and 3CLpro of SARS-CoV-2, respectively. The data suggest that HC exerts collective therapeutic effects against pneumonia caused by SARS-CoV-2 and provides a theoretical basis for further study of the active drug-like ingredients and mechanism of HC in treating pneumonia.


Introduction
Houttuynia cordata Thumb (H. cordata; Saururaceae) is a traditional herbal medicine used in Asian countries, including China, Japan, and Thailand. It exhibits promising antiviral activities towards clinically enveloped viruses, such as influenza virus, herpes simplex virus-1(HSV-1), and human immunodeficiency virus-1 (HIV-1) in vitro [1]. As a time-honored traditional Chinese medicine (TCM), HC has demonstrated a broad range of pharmacological activities for the treatment of inflammatory diseases, especially pulmonary symptoms, i.e., phlegm, dyspnea, lung abscess, and cough, and such infectious diseases as severe acute respiratory syndrome (SARS) [2]. This herb has traditionally been used as one of six principal ingredients in an herbal formula purported to have a preventive effect on SARS-CoV-1 infection [3,4]. As such, it has been recommended to the general public as a

Bioactive Compound Screening and Pharmacokinetic Prediction
HC was searched to identify its active ingredients and pharmacological targets using an in silico approach. As a result, the efficacy of active compounds of HC prediction was established in public databases, including the Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP), the TCM Integrated Database (TCMID), The Encyclopaedia of Traditional Chinese Medicine (ETCM), and the Bioinformatics Analysis Tool for Molecular Mechanism of TCM (BATMAN-TCM) [21]. In addition, absorption, distribution, metabolism, and excretion (ADME) was also employed as a computational evaluation model in pharmacokinetic research to select drug-like compounds [22]. This model consists of criteria such as drug-likeness (DL) and oral bioavailability (OB). These two indices were applied to ascertain whether the compounds have drug-like properties as therapeutic agents and are chemically suitable for drug development [23]. Out of 50 compounds shown in these databases, compounds with DL > 0. 18 and OB > 30% were selected. In total, nine bioactive components were ultimately included in this study and used for the subsequent prediction of compound-related targets.

Potential Targets of HC Active Components
The active components of a drug interact with respective targets to inhibit their biological function. An HC target gene set was acquired by searching several databases: (1) acquiring gene symbols and related information about HC targets from TCMSP; (2) importing selected candidate components into the PubChem database (https://pubchem.ncbi.nlm.nih.gov/ (accessed on 1 August 2021)) to identify relevant targets, and (3) using the ETCM to acquire the target genes associated with the selected active compounds with a score >0.8. The target set was derived after combining the search results and removing duplication and certain bioactive components with suitable targets with a score >0.8 [24].

Identification of Pneumonia-Related Targets Database
Genes related to pneumonia caused by SARS-CoV-2 were screened, selected, and obtained from the Gene Expression Omnibus (GEO), GeneCards database, and other references. First, the dataset was mainly derived from reprocessing high-throughput sequencing data downloaded from the GEO database. The expression profile of GSE152075 involved 430 SARS-CoV-2 infection samples and 54 negative control samples, which were analyzed on the GPL18573 Illumina NexSeq500 platform (Homo sapiens) [25]. The limma package in R Bioconductor was used to identify differentially expressed genes (DEGs) between SARS-CoV-2 infection and negative control, in which the adjusted p-value and [logFC] were determined [26]. The selected criteria of DEGs were thereafter set as p < 0.05 and [logFC] > 1.08 for up-regulated genes and [logFC] < 1.651 for downregulated genes [25]. Second, the targets related to viral pneumonia were obtained by using "Viral pneumonia" as the keyword in the GeneCards database search (https://www.genecards.org/ (accessed on 10 August 2021)). Finally, an intersection between genes retrieved through GeneCards and DEGs of reprocessed data series GSE152075 was obtained as pneumonia-related targets [27].

Construction of PPI Network and Herb-Metabolites-Targets-Disease (HMTD) Network
The HC targets intersecting with pneumonia-related targets were taken as HC-pneumonia common targets, which can be visualized with Venn 2.0. The common targets were imported into the STRING platform (version 11.0). The species was then set to Homo sapience and the minimum required interaction score to the highest confidence of 0.9 in order to retrieve the concise protein-protein interaction (PPI) information for the next step of the analysis. The PPI network was also visualized with Cytoscape software [28].
The Herb-Metabolites-Targets-Disease (HMTD) network was built on the interactions among drug (HC), ingredients, gene symbols, and disease (pneumonia) and then visualized by Cytoscape software. The nodes' varying shapes represent common pneumonia and HC active ingredient targets. The nodes are linked by edges (lines), indicating interactions Viruses 2022, 14, 1588 5 of 20 between nodes. As such, the topological features of a network can be used to predict the targets, while the candidate hub notes analyzed by Cytoscape's Network analyzer tool can be identified by calculating the two topological features of betweenness and degree. Betweenness is the number of shortest paths through the node (the shortest distance between two nodes) and degree is the sum of the number of edges connected the node [29]. Closeness centrality measures how close a node is to all others in the same network [30]. These network centrality indices have been used to define the network properties of drug targets separately or collectively and judge the importance of nodes. Nodes (targets) with higher ranks were considered to have a more critical role within the network [31]. Collectively, the top 15 nodes were screened as the hub genes in the network with the criteria of closeness centrality >0.4 and degree > 4.

GO and KEGG Analysis
Depending on the hub genes (core targets), the Gene Ontology (GO) biological processes, and the Kyoto Encyclopaedia of Genes and Genomes (KEGG), metabolic pathway enrichment analyses were carried out on the pneumonia-HC common targets. GO source and GO enrichment can divide the functions and products of genes into three categories, namely, molecular function (MF), biological process (BP), and cellular components [32]. The biological processes and pathways selected from the analysis of Metascape were colored by cluster ID with the best p-values from each of 20 clusters, wherein nodes that share the same cluster ID are typically close to each other. Enrichment analysis was also carried out in PaGenBase to demonstrate disease targets-organs location. In addition, a further enrichment analysis was performed in DisGeNET to identify the relevant diseases. The results of the KEGG enrichment analysis were used to construct a KEGG pathway network to determine the proteins involved in the treatment effects of HC. Based on the STRING results, the gene-pathway of HC against pneumonia was constructed to delineate the various pathways and key targets in order to explore the potential mechanisms underlying the effect of HC on the treatment of pneumonia [28].

Molecular Docking
Molecular docking is a useful tool to predict and design new drugs. As such, the computational validation of ingredients-targets interactions were confirmed by exploring their binding modes via this process. The computational modeling of intermolecular combinational patterns between target proteins and herb ligands can predict the potential binding modes. Depending on the degree of 67 common targets in the PPI network and the important reference, four active ingredients, including afzelin, apigenin, kaempferol, and quercetin, and six targets, i.e., dipeptidyl peptidase-4 (DPP4), neutrophil elastase (ELANE), Heat shock protein (HSP90AA1), IL6, mitogen-activated protein kinase 1 (MAPK1), and Serpin Family E Member 1 (SERPINE1), were selected to simulate the ingredients-targets interactions for verification of molecular docking. Molecular 2D structures for the molecular ligands and active ingredients were downloaded from PubChem databases. The crystal structures of the key target proteins DPP4 (PDB ID: 4N8D), ELANE (PDB ID: 4WVP), HSP90AA1 (PDB ID: 5J2X), IL6 (PDB ID: 4CNI), MAPK1 (PDB ID: 6QAH), and SERPINE1 (PDB ID: 7AQF) were selected from RCSB PDB (https://www.wwpdb.org/ (accessed on 11 August 2021)). The key target proteins were purposefully selected with a resolution smaller than 2, and their crystals were imported into PyMOL 3.0 software [33]. The active site of the protein is centered on the active amino acid site of the original ligand in the crystal structure, which residue information can be obtained from the literature [34][35][36][37][38][39].
Docking was performed by Autodock Vina 1.1.2, and the molecules with the lowest binding energy in the docking conformation were chosen to observe the binding effect by matching with the original ligands and intermolecular interactions, such as hydrophobic interaction, cation-π, hydrogen bond, anion-π, π-π stacking, salt bridge, and metal complexation [33]. The molecular docking patterns were finally visualized via PyMOL 3.0.

Molecular Dynamics Simulation
The molecular dynamics simulations of these complexes were performed using Gromacs 2020.1, in which the charm36-jul2020 force field was chosen. The complex was solved in TIP3P water and immersed in a dodecahedron box extending to at least 1 nm of the solvent on all sides. The system was neutralized by Na + and Cl − , followed by adding 0.15 M NaCl, which made the system close to the physiological state. The system was minimized by using the steepest descent algorithm for 5000 steps and made a maximum force of less than 1000 kJ/mol/nm. Then, it was equilibrated in a constrained NVT (number of particles, volume, temperature) and NPT (number of particles, pressure, temperature) running for 100 ps. The system was well-equilibrated through NVT and NPT equilibration at 300 K and 1 bar. Finally, MD simulations of the complex were carried out for 100 ns. The Verlet cut-off scheme and a Leap-frog integrator with a step size of 2 fs were applied. The final analysis of molecular dynamics included the root mean square deviation (RMSD) of protein and molecule and the interaction energy between the protein and small molecules, which were calculated by GROMACS 2020.1.

Target Prediction and Analysis of HC
While 50 metabolites were shown in the TCMSP by searching the keywords "Houttuyniae Herba", only 7 satisfied the criteria of OB ≥ 30% and DL ≥ 0.18. Other metabolites were obtained by searching ETCM, BATMAN-TCM, and TCMID. A number of metabolites had already been established as the most effective components of HC throughout the relevant literature [40]. These were also included, even though they did not meet the OB and DL criteria. As such, a total of 20 metabolites were acquired. Depending on these metabolites, the targets for a number of active ingredients of HC were identified by target fishing and by integrating the data acquired from TCMSP, PubChem, and ETCM. The targets of each active ingredient derived from the ETCM database were selected via the screening score ≥0.8. Therefore, only 7 active metabolites were left after searching the targets in those databases, including quercetin, quercitrin, kaempferol, acetyl borneol, decanoic acid, afzelin, and apigenin (Table 1), and 463 targets related to the above seven core active metabolites were identified.

Disease Targets Analysis
Bioinformatics analyses on the expression profile microarray data GSE152075 and GSE1739 containing positive SARS-CoV-2 and negative control samples were performed to identify DEGs between SARS-CoV-2 infection and negative control by the limma package in R Bioconductor. This step identified 9685 DEGs from data series GSE152075 and 1791 DEGs from data series GSE1739. Other disease data sources, such as GeneCards and DisGeNET, were combined with the GEO results to remove duplicates, resulting in the identification of 11,027 targets related to SARS-CoV-2. Since SARS-CoV-2 causes not only pneumonia but also multiple organ failure, neutrophilia, and organ and coagulation dysfunction, pneumonia-related targets were acquired by searching GeneCards with the keywords "SARS-CoV-2" and "pneumonia", resulting in 786 pneumonia-related targets. The intersection between SARS-CoV-2-related targets (11,207) and pneumonia-related targets (786) resulted in 739 elite targets related to pneumonia caused by SARS-CoV-2.

Disease Targets Analysis
Bioinformatics analyses on the expression profile microarray data GSE152075 and GSE1739 containing positive SARS-CoV-2 and negative control samples were performed to identify DEGs between SARS-CoV-2 infection and negative control by the limma package in R Bioconductor. This step identified 9685 DEGs from data series GSE152075 and 1791 DEGs from data series GSE1739. Other disease data sources, such as GeneCards and DisGeNET, were combined with the GEO results to remove duplicates, resulting in the identification of 11,027 targets related to SARS-CoV-2. Since SARS-CoV-2 causes not only pneumonia but also multiple organ failure, neutrophilia, and organ and coagulation dysfunction, pneumonia-related targets were acquired by searching GeneCards with the keywords "SARS-CoV-2" and "pneumonia", resulting in 786 pneumonia-related targets. The intersection between SARS-CoV-2-related targets (11,207) and pneumonia-related targets (786) resulted in 739 elite targets related to pneumonia caused by SARS-CoV-2.

Disease Targets Analysis
Bioinformatics analyses on the expression profile microarray data GSE152075 and GSE1739 containing positive SARS-CoV-2 and negative control samples were performed to identify DEGs between SARS-CoV-2 infection and negative control by the limma package in R Bioconductor. This step identified 9685 DEGs from data series GSE152075 and 1791 DEGs from data series GSE1739. Other disease data sources, such as GeneCards and DisGeNET, were combined with the GEO results to remove duplicates, resulting in the identification of 11,027 targets related to SARS-CoV-2. Since SARS-CoV-2 causes not only pneumonia but also multiple organ failure, neutrophilia, and organ and coagulation dysfunction, pneumonia-related targets were acquired by searching GeneCards with the keywords "SARS-CoV-2" and "pneumonia", resulting in 786 pneumonia-related targets. The intersection between SARS-CoV-2-related targets (11,207) and pneumonia-related targets (786) resulted in 739 elite targets related to pneumonia caused by SARS-CoV-2.

Disease Targets Analysis
Bioinformatics analyses on the expression profile microarray data GSE152075 and GSE1739 containing positive SARS-CoV-2 and negative control samples were performed to identify DEGs between SARS-CoV-2 infection and negative control by the limma package in R Bioconductor. This step identified 9685 DEGs from data series GSE152075 and 1791 DEGs from data series GSE1739. Other disease data sources, such as GeneCards and DisGeNET, were combined with the GEO results to remove duplicates, resulting in the identification of 11,027 targets related to SARS-CoV-2. Since SARS-CoV-2 causes not only pneumonia but also multiple organ failure, neutrophilia, and organ and coagulation dysfunction, pneumonia-related targets were acquired by searching GeneCards with the keywords "SARS-CoV-2" and "pneumonia", resulting in 786 pneumonia-related targets. The intersection between SARS-CoV-2-related targets (11,207) and pneumonia-related targets (786) resulted in 739 elite targets related to pneumonia caused by SARS-CoV-2.

Disease Targets Analysis
Bioinformatics analyses on the expression profile microarray data GSE152075 and GSE1739 containing positive SARS-CoV-2 and negative control samples were performed to identify DEGs between SARS-CoV-2 infection and negative control by the limma package in R Bioconductor. This step identified 9685 DEGs from data series GSE152075 and 1791 DEGs from data series GSE1739. Other disease data sources, such as GeneCards and DisGeNET, were combined with the GEO results to remove duplicates, resulting in the identification of 11,027 targets related to SARS-CoV-2. Since SARS-CoV-2 causes not only pneumonia but also multiple organ failure, neutrophilia, and organ and coagulation dysfunction, pneumonia-related targets were acquired by searching GeneCards with the keywords "SARS-CoV-2" and "pneumonia", resulting in 786 pneumonia-related targets. The intersection between SARS-CoV-2-related targets (11,207) and pneumonia-related targets (786) resulted in 739 elite targets related to pneumonia caused by SARS-CoV-2.

Disease Targets Analysis
Bioinformatics analyses on the expression profile microarray data GSE152075 and GSE1739 containing positive SARS-CoV-2 and negative control samples were performed to identify DEGs between SARS-CoV-2 infection and negative control by the limma package in R Bioconductor. This step identified 9685 DEGs from data series GSE152075 and 1791 DEGs from data series GSE1739. Other disease data sources, such as GeneCards and DisGeNET, were combined with the GEO results to remove duplicates, resulting in the identification of 11,027 targets related to SARS-CoV-2. Since SARS-CoV-2 causes not only pneumonia but also multiple organ failure, neutrophilia, and organ and coagulation dysfunction, pneumonia-related targets were acquired by searching GeneCards with the keywords "SARS-CoV-2" and "pneumonia", resulting in 786 pneumonia-related targets. The intersection between SARS-CoV-2-related targets (11,207) and pneumonia-related targets (786) resulted in 739 elite targets related to pneumonia caused by SARS-CoV-2.

Disease Targets Analysis
Bioinformatics analyses on the expression profile microarray data GSE152075 and GSE1739 containing positive SARS-CoV-2 and negative control samples were performed to identify DEGs between SARS-CoV-2 infection and negative control by the limma package in R Bioconductor. This step identified 9685 DEGs from data series GSE152075 and 1791 DEGs from data series GSE1739. Other disease data sources, such as GeneCards and DisGeNET, were combined with the GEO results to remove duplicates, resulting in the identification of 11,027 targets related to SARS-CoV-2. Since SARS-CoV-2 causes not only pneumonia but also multiple organ failure, neutrophilia, and organ and coagulation dysfunction, pneumonia-related targets were acquired by searching GeneCards with the keywords "SARS-CoV-2" and "pneumonia", resulting in 786 pneumonia-related targets. The intersection between SARS-CoV-2-related targets (11,207) and pneumonia-related targets (786) resulted in 739 elite targets related to pneumonia caused by SARS-CoV-2.

Herb-Ingredients-Targets-Disease Network of HC Analysis
The intersection between HC-related targets (463) and pneumonia-related elite targets (739) resulted in an HC-pneumonia common target set with 67 genes. This common target set was imported into Cytoscape v3.5.0 to construct an Herb-Metabolites-Targets-Disease (HMTD) network, as shown in Figure 2. This consisted of six metabolites assigned to 67 targets, indicating HMTD interactions. The node can be designed as a hub node if the degree, betweenness, and closeness satisfy specific criteria, such as the median of the corresponding parameters. The screening of important metabolites and core targets was carried out based on the criteria of SUID > 70, Closeness Centrality > 0.4 and Degree > 4, resulting in four ingredients (apigenin, quercetin, afzelin, and kaempferol) and 21 core common targets (Table 2). Consequently, these important metabolites might be crucial active compounds of HC targeting 21 genes, which could be verified through molecular docking or further experiments. The HC-pneumonia common target set was imported into STRING to remove unconnected targets, and a PPI network with a confidence score set to 0.9 or higher was gained. PPI information from the STRING platform was input to Cytoscape software to construct a PPI network based on the common targets shown in Figure 3. The size of target nodes was consistent with the degree, and the nodes with pink color were deemed to be important targets [41]. According to the degree and combined score, the top 21 common targets shown in Table 2 are involved in the effects of HC treatment on pneumonia.

GO and KEGG Enrichment Analysis
According to the results of GO enrichment analyses in Metascape, the genes were enriched in different GO terms, and the top 20 GO terms in the three categories were selected to construct connections within the signaling network ( Figure 4). The network was colored by cluster ID with the best p-values from each of 20 clusters, wherein nodes that share the same cluster-ID are typically close. The top signaling pathways mainly include interleukin-4 and interleukin-13 signaling, the AGE-RAGE signaling pathway in diabetic complication, positive regulation of cell migration, spinal cord injury, IL17 signaling pathway, reactive oxygen species metabolic process, leukocyte activation involved in immune response, Th17 cell differentiation, cytokines and inflammatory response, and epithelial cell migration ( Figure 4). The x-axis in the bar chart represents log 10 (p-value) and the y-axis the GO term. The enrichment analysis in PaGenBase demonstrated disease targets-organs location, such as lung, smooth muscle, cardiac myocytes, bone marrow, bronchial epithelial cells, liver cells, and spleen. The analyses of the most-associated diseases showed immunosuppression, fatty acid disease, respiratory distress syndrome, pneumonitis, endothelial dysfunction, respiratory syncytial virus infection, liver failure, middle cerebral artery occlusion, bacterial infections, acute myocardial infarction, lung diseases, cardiac arrest, myocardial ischemia, and herpes simplex infections. To analyze the significance and importance of key targets in the pathways involved in the treatment effect of pneumonia, 10 key pathways, determined according to gene counts and adjusted p values from the KEGG enrichment analysis and related targets, were used to construct a KEGG key pathway network ( Figure 5). The construction of gene-KEGG key pathways demonstrated that the targets MAPK1, MAPK3, IL6, PIK3CA, AKT1, EGFR, TNF, and STAT3 involved more than four signaling pathways, including PI3K-Akt, Jak-STAT, MAPK, and NF-kB (Table 3). As such, HC could target multiple functional and biological factors in pneumonia. However, the effects and profound influence required further validation [28].   Table 3.

Docking and Molecular Dynamics Simulation Analysis of Ingredients-Targets
Considering the integration of the results from the PPI network and HMTD network, the key targets in the pathways mentioned above and the nodes with high degrees represent the key targets. Therefore, molecular docking validation was performed based on the pneumonia-related targets and selected ingredients from the HMTD network. The selected active ingredients included quercetin (CAS no. 117-39-5), kaempferol (CAS no. 520-18-3), afzelin (CAS no. 482-39-3), and apigenin (CAS no. 520- . The protein structures of key targets were acquired online from RCSB PDB, including DPP4, ELANE, HSP90AA1, IL6, MARK1, and SERPINE1, based on STRING interaction analysis and importance reference. Docking analysis of the metabolites and proteins above showed the docking patterns and binding affinities ( Table 4). The docking results were represented on the molecular surface to reflect the topical details of the binding sites. The residues were marked on the protein surface, and hydrogen bonds were shown as solid lines ( Figure 6). The binding affinities of all docking patterns were less than −6 kcal/mol, indicating a stable binding between active ingredients and protein targets. The binding affinities are listed in Table 4, and the binding configuration is shown in Figure 6. The affinity energy of the best mode, apigenin, is −9.4 kcal/mol.   Table 4.
In order to verify the stability of the docking structures, we selected DPP4-kaempferol, MAPK1-afzelin, SERPINE1-apigenin and SERPINE1-quercetin complexes for dynamic simulation analysis. As shown in Figure 7, RMSD of proteins and small molecules in the complex structures remained relatively stable during the simulation, especially MAPK1-afzelin, SERPINE1-apigenin and SERPINE1-quercetin complexes. The RMSD of kaempferol in DPP4-kaempferol varied greatly, and it was stable at 1 nm com-  Table 4.
In order to verify the stability of the docking structures, we selected DPP4-kaempferol, MAPK1-afzelin, SERPINE1-apigenin and SERPINE1-quercetin complexes for dynamic simulation analysis. As shown in Figure 7, RMSD of proteins and small molecules in the complex structures remained relatively stable during the simulation, especially MAPK1-afzelin, SERPINE1-apigenin and SERPINE1-quercetin complexes. The RMSD of kaempferol in DPP4-kaempferol varied greatly, and it was stable at 1 nm compared with the initial docking structure, indicating that the position of kaempferol changed significantly during the simulation process. This position change occurred rapidly and remained stable at the new binding position. The average interaction energy of DPP4-kaempferol, MAPK1-afzelin, SERPINE1-apigenin and SERPINE1-quercetin complexes was −107.34 kJ/mol, −180.39 kJ/mol, −175.82 kJ/mol, and −183.77 kJ/mol. In conclusion, we used docking and dynamics simulation to explore ingredient-target prediction. Our results showed that the DPP4-kaempferol, MAPK1-afzelin, SERPINE1apigenin, and SERPINE1-quercetin complexes had good docking fractions, that the protein and small-molecule positions were stable during the simulation process, and that the interaction energy was lower than −100 kJ/mol. Accordingly, our studies provide a reference for subsequent experimental design.

Discussion
While previous studies demonstrated that HC exhibited antiviral activity through inhibiting SARS-CoV-1 3CLPRO and RdRp and also stimulated the proliferation of CD4 + and CD8 + T cells in vitro [3], the ability of HC to inhibit SARS-CoV-2 infection and the similarities between the inhibitory mechanisms of SARS-CoV-2 and SARS-CoV-1 remain unknown. Although several clinical studies have confirmed that prescriptions or formulas (Lian Hua Qing Wen capsule) containing HC are effective for the treatment of COVID-19 [42], further research is necessary to investigate the mechanism of action of HC and apply the pharmacology of TCM network for the predictive analysis. Thus, this study aimed to explore the pharmacological mechanism of HC on pneumonia caused by SARS-CoV-2 using network pharmacology and molecular docking.
TCM considers an individual or patient an integrative complex with dynamic states, demonstrating multiple biological targets and focusing on integral therapeutic efficacies [43]. Inflammation has been the pathophysiological mechanism behind many chronic diseases, including cytokines, nitric oxide (NO), lipid mediators, G prostaglandins, and leukotrienes produced by macrophages, neutrophils, and other inflammatory cells [44,45]. This study included seven core metabolites, following diligent screening and searching in the references, including quercetin, quercitrin, kaempferol, acetyl borneol, decanoic acid, afzelin, and apigenin. These flavonoids are known to be large entities of plant constituents and possess anti-inflammatory activity [15]. Some flavonoids have been shown to attenuate lung inflammatory response strongly. For example, quercetin was previously reported to attenuate lipopolysaccharide (LPS)-induced lung inflammation in mice by oral administration [14], while afzelin isolated from methanol extract of HC was demonstrated to regulate both mitophagy and mitochondrial biogenesis through Rev-Erb-/phosphor-AMPK/SIRT1 signaling [40]. In addition, afzelin can inhibit mitochondrial dysfunction induced by excessive oxidative stress and attenuate the reduction of mitochondrial GDH activity and hepatic ATP production in LPS-induced hepatic injury [40]. Kaempferol has been proven to protect against H9N2 swine influenza virus infection and can ameliorate virus-induced acute lung injury by inactivation of TLR4/MyD88-mediated NF-kB and MAPK signaling pathways [46]. Kaempferol can also inhibit the release of TNF-α, IL1β, IL6, and IL18 and suppress the activation of NF-kB and AKT, thus attenuating cardiac fibroblast inflammation [47]. Apigenin, a plant-derived flavonoid, possesses anti-carcinogenic, antioxidant, anti-inflammatory, and anti-mutagenic properties [48]. Apigenin can react to the Nrf2 gene, which encodes a key transcription factor to regulate the antioxidative defense system, and was also a potent inhibitor of SARS-CoV 3CLpro [47]. The targets of these active ingredients were collected from a different database.
Depending on the above active metabolites, target searching in the databases of TCMSP, PubChem, and ETCM resulted in the acquisition of 463 targets related to the above seven core active metabolites. The disease targets were obtained by identifying key genes in SARS-CoV-2 infection and uncovering their potential functions by re-processing the expression profiling of high-throughput sequencing of GSE152075 from the GEO database to form a solid basis for the ensuing analysis. The bioinformatics analyses generated 11,207 DEGs, contributing to our understanding of the molecular mechanism underlying the advancement of SARS-CoV-2 infection.
Since SARS-CoV-2 causes multiple organ failure [49], the intersection between DEGS acquired from the GEO series and pneumonia-related genes from GeneCards and Dis-GeNET isolated pneumonia caused by SARS-CoV-2, resulting in 739 elite targets. These were matched and mapped to obtain 67 common HC-pneumonia targets. A PPI network of the common targets and screened nodes according to observed gene count > 50 resulted in 21 core target networks to verify targets associated with HC ingredients. According to the combined score of each node, the top 25 nodes were mainly NOS3, MTOR, SERPINE1, PIK3CA, HSP90AA1, STAT3, INS, IL6, IL-1β, TNF, IL10, VEGFA, CDK1, and GSK3B. Most of these genes are associated with inflammation, hypoxia, and angiogenesis [50]. This result aligned with the findings of previous studies of HC extract and confirmed that HC can down-regulate specific inflammatory mediators, such as TNF-α, IL6, prostaglandin E2 (PGE2), and nitric oxide (NO) production in the cells, inducible nitric oxide synthase (iNOS), and cyclooxygenase-2 (COX-2) expression [7].
The pathogenesis of inflammatory diseases is associated with the overproduction of the above mediators. NO by endothelial nitric oxide synthase (NOS3) is implicated in vascular smooth muscle relaxation and mediates vascular endothelial growth factor (VEGF)-induced angiogenesis in coronary vessels, which may explain the many complications occurring in COVID-19 patients [51]. The mammalian target of rapamycin (mTOR) entails the downstream of PI3K-Akt to regulate cell growth and proliferation, cell survival, protein synthesis, and transcription [52]. HC may inhibit PI3K/Akt/mTOR and ERK1/2 signaling pathways in human lung cancer cells [53], while IL17A can inhibit PI3K/Akt/mTOR-mediated autophagy, which causes lung inflammation and fibrosis [54]. However, alternative evidence suggests that activation of the PI3K/Akt/mTOR pathway may contribute to pulmonary fibrosis and lung injury by regulating lung fibroblasts and epithelial cells [55]. For STAT3, the activation of JAK-STAT signaling can lead to fibrosis in many organs, i.e., the lung [47]. STAT3 is a potential molecular target for clinical syndromes characterized by systematic inflammation in COVID-19 in a large-scale transcriptional study [56]. Therefore, GO and KEGG enrichment analysis of the core targets was applied to explore the underlying mechanism of HC.
The GO enrichment results confirmed that HC mainly regulates biological processes in response to stimulus, regulation of the cellular process, response to chemicals, response to stress, cell communication, signal transduction regulation of macromolecule metabolic process, and regulation of nitrogen compound metabolic processes. The results of KEGG enrichment uncovered the numerous signaling pathways involved in the development and progression of pneumonia, including PI3K-Akt, HIF-1, IL-17, TNF, TLR, JAK-STAT, NOD-like receptor, or MAPK. Previous studies further suggest that the anti-inflammatory properties of HC extract may arise from the inhibition of pro-inflammatory mediators by suppression of NF-κB, and MAPK signaling pathways by binding their key proteins in pathways with HC active metabolites, as outlined above [9,39,41].
Based on the PPI, HMTD, and KEGG enrichment analysis, 6 of 21 core targets, namely, DPP4, ELANE, HSP90AA1, IL6, MAPK1, and SERPINE1, were selected to dock with four main ingredients identified in a series of analyses of this study. As such, the docking results agreed with the intermolecular interactions. Despite extensive prior efforts to elucidate the metabolism and effects of oral administration of HC in patients in vivo, the mechanism remained largely unclear. DPP4, HSP90AA1, and SERPINE1 are related to immune response and linked to macrophages which promote the generation of inflammatory factors, such as TNF-α and IL6 [57]. DPP4 is a member of serine peptidases known as adenosine deaminase complexing protein 2 clusters of differentiation 26 (CD26) associated with immune regulation, apoptosis, and signal transduction. However, the main receptor of SARS-CoV is angiotensin-converting enzyme 2 (ACE2) or CD209L, whereas MERS-CoV uses DPP4 (also known as CD26) as the major receptor [58]. A correlation between DPP4 and ACE2 found that both membrane proteins can facilitate virus entry. Therefore, DPP4 was speculated to be a co-receptor to facilitate SARS-CoV-2 infection since DPP4 can be found in lung cells [59]. The co-receptors of ACE2 and DPP4 to the spike glycoprotein postulate that different human coronaviruses target similar cell types across different human organs [60]. Heat shock proteins (HSPs), also known as stress proteins, are divided into the HSP70 family, HSP90 family, HSP 100 family, and so on. HSP90AA1 is a molecular chaperone protein with inhibitors that can induce hepatic stellate cell apoptosis through neurophospholipase or NK-κB, depending on the mechanism. A previous study found that Jin Hua Qing Gan Granules regulate multiple signaling pathways via binding targets, such as PTGS2, HSP90AA1, and NCOA2, to prevent COVID-19 [61].
ELANE genes can express up to a five-fold level against SARS-CoV-1 infection, causing lung proinflammatory cytokines [62]. Recent research found strong evidence that the ELANE gene mainly enhances proinflammatory cytokines and can subsequently cause epithelial cell injuries in cystic fibrosis patients [62,63]. The suppression of the ELANE gene can therefore reduce the production of proinflammatory cytokines, resulting in improved pulmonary function [63]. As such, the inhibition of ELANE directly protects the lung and reduces lung inflammatory cell infiltration to improve the success rate of coronavirus patients [62]. MAPK1 is the critical target participating in the MAPK signaling and PI3K-Akt signaling pathways. The inhibition of MAPK1 can result in inhibiting the above two pathways in LPS-induced ALI in mice [64]. SERPINE1 is a serpin peptidase inhibitor whose increased expression resulted in a lower survival rate. Furthermore, increased expression of SERPINE1 was associated with the activation of the PI3K-Akt pathway [16]. The results from PPI and HMTD indicated that the six intersection genes above could prove to be potent pharmacological targets of HC against COVID-19.
Molecular docking of 6 out of 21 core targets and four screened active ingredients was applied to validate the results of the network pharmacology. The active ingredients which target these proteins include quercetin, kaempferol, afzelin, and apigenin. Docking within the docking pockets between active ingredients and target proteins was visualized by AutoDock Vina and PyMOL software. The binding affinities of the docking results ranged from −6.7 to −9.4 kcal/mol, indicating stable binding. Afzelin showed the greatest binding affinities with MAPK1 (−9.4 kcal/mol), followed by DPP4 (−8.6 kcal/mol). In contrast, MAPK showed a better average binding affinity with these four active ingredients, indicating that HC protected against acute lung injury, mainly through the suppression of the MAPK/NF-kB pathway [5]. The docking pose of apigenin shows H-bonds between the aromatic region and residues PHE-129, ILE-133, ASN-82, ASP-106, ILE-84, ASN-158, and THR-190, establishing a stacking interaction with GLN-132 [65]. However, the docking results reflect possible treatment mechanisms and guide herb-disease validation via cells and animal experiments.
Network pharmacology has been widely used for TCM mechanism research owing to its efficacy in analyzing the complicated relationships among multiple ingredients and multiple disease targets. Docking technology can visualized the binding modes of active ingredients with disease-related key target proteins and provide guidance for researchers' selections of active ingredients for in vivo or in vitro experiments [66]. As shown in the binding affinity result, multiple ingredients can bind the same protein, which may result in synergistic effects. This represents a significant challenge for current network pharmacology. In addition, further animal and cell models are needed to verify the relevant pathways and targets.
This study confirmed the potent therapeutic effect of HC, a time-honored herb widely used in Asian countries to treat pneumonia. The potential mechanisms of HC were revealed by employing both network pharmacology and molecular computational analyses. Our results offer a very different perspective in terms of modern pharmacological mechanisms which may assist in the global fight against the COVID-19 pandemic. However, while HC is an effective herb for pneumonia treatment, the optimal dose for inducing remission with low toxicity needs to be determined. Moreover, further animal and cell models are necessary to verify the relevant pathways and targets.
Author Contributions: J.L., S.J., H.S. and J.W. participated in the design of this study. J.L. and Y.Y. performed the statistical analysis. J.L. and H.S. drafted the manuscript, while G.S., S.J. and H.S. revised the manuscript. S.Y. performed the molecular dynamics simulation. All data were generated in-house, and no paper mill was used. All authors have read and agreed to the published version of the manuscript.
Funding: This study is supported by the grants from the H2020 Marie Skłodowska-Curie Actions and Enterprise Ireland (No. 713654), the European Union's Horizon 2020 Research and Innovation Programme under Marie Skłodowska-Curie grant agreement (No. 754489), the Science Foundation Ireland grant (No. 13/RC/2094) and the National Natural Science Foundation of China (No. 82161138002). The funding sources had no role in the study design, data collection, analysis, interpretation, or writing of the report.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The raw data of this paper are available on request.